Data preparation

Data are cleaned and saved in RData format for the subsequent analyses.

# Step 1: load the original data
con = dbConnect(SQLite(),
                dbname="VariantSCN5A-second-revision.db")
alltables = dbListTables(con)
my.data <- dbReadTable(con, 'VariantSCN5A')
my.data[my.data=='NA'] <- NA
data<-my.data
dbDisconnect(con)

dim(data) # 2417 by 29
names(data)

# clean `resnum`
table(data$resnum)
sum(is.na(data$resnum)) # 51 * + 1246 missing = 1297 abnormal
sum(data$resnum=="") # 1246 missing

# convert strings to numeric values
data$resnum<-suppressWarnings(as.integer(data$resnum))
sum(is.na(data$resnum)) # 1300 missing, three more 

# replace missing resnum with the number after the first capital letters?
for (variant in data$var){
  if (is.na(data[data$var==variant, "resnum"])){
    suppressWarnings(data[data$var==variant, "resnum"]<-as.integer(strsplit(variant,"[A-Z]")[[1]][2]))
  }
}

sum(is.na(data$resnum)) # 45 unsuccessful conversion

data[is.na(data$resnum),c("var","resnum")]

# manually get resnum from var among variants of interest
data[data$var=="E444fsX14","resnum"]<-444
data[data$var=="F1775fsX1785","resnum"]<-1775
data[data$var=="L1579fsX53","resnum"]<-1579
data[data$var=="p.1795insD","resnum"]<-1795
data[data$var=="p.G1031fsX27","resnum"]<-1031
data[data$var=="p.G1748del","resnum"]<-1748
data[data$var=="p.G1758del","resnum"]<-1758
data[data$var=="p.L1339del","resnum"]<-1339
data[data$var=="p.M741_T742delinsI ","resnum"]<-741
data[data$var=="Y1795insD","resnum"]<-1795
data[data$var=="p.A586_L587del","resnum"]<-586
data[data$var=="p.E1939_E1943del","resnum"]<-1939
data[data$var=="p.E1072del","resnum"]<-1072
data[data$var=="p.E1064del","resnum"]<-1064
data[data$var=="p.F1465_L1480dup","resnum"]<-1465
data[data$var=="p.M1875dup","resnum"]<-1875
data[data$var=="p.N507_L515dup","resnum"]<-507
data[data$var=="p.N291_S293dup","resnum"]<-291
data[data$var=="p.Pro2006del","resnum"]<-2006
data[data$var=="p.Gln1000del","resnum"]<-1000
data[data$var=="p.S1970_S1972del","resnum"]<-1970
data[data$var=="p.T290_G292del","resnum"]<-290
data[data$var=="p.I1758del","resnum"]<-1758


data[is.na(data$resnum),c("var","resnum")]

# after the conversion, how many missing in resnum?
sum(is.na(data$resnum)) # 23

# are resnum unique? No
length(unique(data$resnum))

# are nativeAA and mutAA different when resnum is the same?
duplicates <- data[duplicated(data$resnum)|duplicated(data$resnum,fromLast = T),c("var","resnum","nativeAA","mutAA")]

# it appears so  
head(duplicates[order(duplicates$resnum),])

# nativeAA
table(data$nativeAA)
sum(is.na(data$nativeAA)) # 168 missing

# mutAA
table(data$mutAA)
sum(is.na(data$mutAA)) # 168 missing

# create a composite index for merging
data$res_nm <- paste0(as.character(data$resnum),data$nativeAA,data$mutAA)

# Step 2: load five supplementary data
# load data
bcl  <- read.csv("covariates/scn5a_bcl_features.txt", header = TRUE, sep = "\t")
sasa <- read.csv("covariates/5a_5x0m_sasa.csv", header = TRUE)
pph2<-read.csv("covariates/pph2-20180329.csv", strip.white = TRUE)
prov<-read.csv("covariates/provean_revised.csv", strip.white = TRUE)
pam <- read.csv("covariates/sin_pam_energy_funcdist.csv", header = TRUE)

# Step 3: examine supplementary data
# pam
# pam has the same covariates as in the original data
sum(names(pam)%in%names(data)) == length(names(pam))
# seven times duplicated observations
nrow(pam) / nrow(unique(pam))
pam.u <- unique(pam)
head(pam.u)
data[data$var == "A2T","pamscore"]
# no need to merge with pam

# pph2 
dim(pph2)
# all unique
dim(unique(pph2))
# contains four more covariates but no `var`
sum(names(pph2)%in%names(data)) == length(names(pph2))
names(pph2)%in%names(data)
names(pph2)
head(pph2)
# each residue number has 19 rows
length(unique(pph2$resnum))
pph2[pph2$resnum==1,]
table(data$resnum)
# create a composite index for merging
pph2$res_nm <- paste0(as.character(pph2$resnum),pph2$nativeAA,pph2$mutAA)

# sasa
dim(sasa)
# all unique
dim(unique(sasa))
head(sasa)
# unique resnum
length(unique(sasa$resnum)) == nrow(sasa)

# bcl
dim(bcl)
# eight duplicates
dim(unique(bcl))
head(bcl)
bcl[duplicated(bcl),]
bcl <- distinct(bcl)
dim(bcl)
# no duplicates
dim(unique(bcl))
head(bcl)

# prov
dim(prov)
dim(unique(prov))
# 15 duplicates
prov <- distinct(prov)
dim(unique(prov))[1] == nrow(prov)
head(prov)
# create a composite index for merging
prov$res_nm <- paste0(as.character(prov$resnum),prov$nativeAA,prov$mutAA)


# Step 4: merge data with 4 additional data sets
# pph
data.pph <- merge(data,pph2,by="res_nm",all.x=T)
dim(data.pph)
# drop the resnum, mutAA, nativeAA from pph2
data.pph<-data.pph[,-which(names(data.pph) %in% c("resnum.y","mutAA.y","nativeAA.y"))]

# rename resnum
colnames(data.pph)[colnames(data.pph)=="resnum.x"] <- "resnum"
colnames(data.pph)[colnames(data.pph)=="mutAA.x"] <- "mutAA"
colnames(data.pph)[colnames(data.pph)=="nativeAA.x"] <- "nativeAA"
names(data.pph)

# sasa
data.pph.sasa <- merge(data.pph,sasa,by="resnum",all.x=T)
dim(data.pph.sasa)
colnames(data.pph.sasa)

# bcl
data.pph.sasa.bcl <- merge(data.pph.sasa,bcl,by="var",all.x=T)
dim(data.pph.sasa.bcl)
names(data.pph.sasa.bcl)

# prov
data.pph.sasa.bcl.prov <- merge(data.pph.sasa.bcl,prov,by="res_nm",all.x = T)
dim(data.pph.sasa.bcl.prov)
names(data.pph.sasa.bcl.prov)

# drop the resnum, mutAA, nativeAA from prov
data.pph.sasa.bcl.prov<-data.pph.sasa.bcl.prov[,-which(names(data.pph.sasa.bcl.prov) %in% c("resnum.y","mutAA.y","nativeAA.y"))]

# rename resnum
colnames(data.pph.sasa.bcl.prov)[colnames(data.pph.sasa.bcl.prov)=="resnum.x"] <- "resnum"
colnames(data.pph.sasa.bcl.prov)[colnames(data.pph.sasa.bcl.prov)=="mutAA.x"] <- "mutAA"
colnames(data.pph.sasa.bcl.prov)[colnames(data.pph.sasa.bcl.prov)=="nativeAA.x"] <- "nativeAA"
names(data.pph.sasa.bcl.prov)
dim(data.pph.sasa.bcl.prov)
head(data.pph.sasa.bcl.prov)
mydata <- data.pph.sasa.bcl.prov

# Step 5: exclusion criteria
# drop zero total carriers

mydata$gnomAD<-as.integer(mydata$gnomAD)
sum(is.na(mydata$gnomAD))
mydata$gnomAD[is.na(mydata$gnomAD)] <- 0

mydata$total_carriers<-mydata$unaff+mydata$lqt3+mydata$brs1+mydata$gnomAD
mydata.total<-mydata[mydata$total_carriers>0, ]
dim(mydata.total)
names(mydata.total)

# keep only missense, aadel and aains mututation types
mydata.total <- mydata.total[mydata.total$mut_type %in% c("missense","aadel", "aains"),]
dim(mydata.total)

# rename the final data set without clinical covariates
d <- mydata.total
d<-d[,!(names(d) %in% "aaneigh_rank_inc_loops")]
d<-d[!is.na(d$resnum),]
dim(d)
names(d)
save(d,file="BrS1_data.RData")

EM algorithm

Step 3 Get posterior \(\alpha_{1i}\) and \(\beta_{1i}\)

Variant-specific Beta-Binomial posterior \(\alpha_{1i}\) and \(\beta_{1i}\) based on Beta prior

\[\alpha_{1i}=\alpha_0+x_i,\beta_{1i}=\beta_0+n_i-x_i \]

Back to top

Step 4 Get posterior mean \(\tilde\mu_i\) and variance \(\tilde\sigma^2_i\)

The mean \(\tilde\mu_i\) of posterior distribution is

\[\tilde \mu_i=\frac{\alpha_{1i}}{\alpha_{1i}+\beta_{1i}} \]

The variance \(\tilde\sigma^2_i\) of posterior distribution is

\[\tilde\sigma^2_i=\frac{\alpha_{1i}\beta_{1i}}{(\alpha_{1i}+\beta_{1i})^2(\alpha_{1i}+\beta_{1i}+1)} \]

which is raised at 20th percentile.

Back to top

Step 5 Calculate penetrance density from \(\tilde\mu_i\) using function funcdist

Note that only 774 out of 1439 (53.8%) variants fall within structured regions of the protein.

## [1] "p.Gln1000del"
## [1] "Q1000L"
## [1] "P1002S"
## [1] "C1004R"
## [1] "I1005T"
## [1] "I1005V"
## [1] "A1006S"
## [1] "T1007I"
## [1] "T1007N"
## [1] "P1008S"
## [1] "P1011L"
## [1] "P1011S"
## [1] "P1014S"
## [1] "E1015K"
## [1] "T1016M"
## [1] "K1018E"
## [1] "T101I"
## [1] "P1021S"
## [1] "R1023C"
## [1] "R1023H"
## [1] "R1023P"
## [1] "K1024R"
## [1] "E1025A"
## [1] "R1027P"
## [1] "R1027Q"
## [1] "R1027W"
## [1] "E1029K"
## [1] "E1032D"
## [1] "E1032K"
## [1] "Q1033R"
## [1] "P1034T"
## [1] "G1035V"
## [1] "G1040R"
## [1] "D1041G"
## [1] "D1041N"
## [1] "E1043K"
## [1] "V1045M"
## [1] "P1048S"
## [1] "R104G"
## [1] "R104Q"
## [1] "R104W"
## [1] "A1050T"
## [1] "V1051A"
## [1] "A1052D"
## [1] "E1053K"
## [1] "D1055G"
## [1] "T1056A"
## [1] "E1061D"
## [1] "D1062H"
## [1] "E1063G"
## [1] "p.E1064del"
## [1] "S1066G"
## [1] "L1067R"
## [1] "G1068A"
## [1] "G1068D"
## [1] "T1069M"
## [1] "S106T"
## [1] "E1071K"
## [1] "p.E1072del"
## [1] "S1074G"
## [1] "S1074R"
## [1] "S1079F"
## [1] "S1079T"
## [1] "S1079Y"
## [1] "V1082A"
## [1] "S1083C"
## [1] "G1084D"
## [1] "G1084R"
## [1] "G1084S"
## [1] "A1088T"
## [1] "A1088V"
## [1] "P1090L"
## [1] "P1090Q"
## [1] "D1091A"
## [1] "D1091Y"
## [1] "W1095C"
## [1] "S1096C"
## [1] "S1096G"
## [1] "Q1097H"
## [1] "V1098L"
## [1] "V1098M"
## [1] "N109K"
## [1] "T10S"
## [1] "A1100T"
## [1] "A1100V"
## [1] "A1102T"
## [1] "S1103F"
## [1] "S1103Y"
## [1] "E1105V"
## [1] "A1106T"
## [1] "E1107K"
## [1] "S1109G"
## [1] "A110T"
## [1] "A1113T"
## [1] "A1113V"
## [1] "D1114E"
## [1] "D1114N"
## [1] "W1115R"
## [1] "R1116Q"
## [1] "R1116W"
## [1] "A1121V"
## [1] "A1125G"
## [1] "A1125T"
## [1] "A1125V"
## [1] "G1129S"
## [1] "Y112C"
## [1] "E1130K"
## [1] "T1131I"
## [1] "T1131S"
## [1] "P1132S"
## [1] "D1134E"
## [1] "D1134N"
## [1] "S1135I"
## [1] "C1136Y"
## [1] "V113A"
## [1] "V113I"
## [1] "S1140T"
## [1] "T1147I"
## [1] "T1147N"
## [1] "T1147S"
## [1] "A1148S"
## [1] "A1148T"
## [1] "Q1153H"
## [1] "I1154N"
## [1] "P1155S"
## [1] "D1156G"
## [1] "G1158S"
## [1] "S115G"
## [1] "D1163E"
## [1] "D1163G"
## [1] "D1163Y"
## [1] "P1164T"
## [1] "E1165D"
## [1] "E1165Q"
## [1] "D1166N"
## [1] "C1167Y"
## [1] "F1168L"
## [1] "T1169I"
## [1] "V1173D"
## [1] "R1174G"
## [1] "R1174W"
## [1] "R1175H"
## [1] "P1177L"
## [1] "C1178Y"
## [1] "A1180V"
## [1] "V1181A"
## [1] "V1181L"
## [1] "V1181M"
## [1] "T1183I"
## [1] "A1186T"
## [1] "P1187Q"
## [1] "K1189T"
## [1] "V1190F"
## [1] "R1193Q"
## [1] "R1193W"
## [1] "L1194M"
## [1] "R1195H"
## [1] "R1195S"
## [1] "Y1199S"
## [1] "P119L"
## [1] "P119S"
## [1] "S11R"
## [1] "H1200R"
## [1] "H1200Y"
## [1] "I1201M"
## [1] "V1202M"
## [1] "E1208K"
## [1] "T1209R"
## [1] "F1210S"
## [1] "p.I1212del"
## [1] "M1214T"
## [1] "I1215V"
## [1] "L1216V"
## [1] "S1218I"
## [1] "S1218T"
## [1] "S1219N"
## [1] "R121Q"
## [1] "R121W"
## [1] "G1220E"
## [1] "A1221V"
## [1] "L1222R"
## [1] "E1225K"
## [1] "G1225K"
## [1] "Y1228C"
## [1] "Y1228F"
## [1] "Y1228H"
## [1] "E1230K"
## [1] "E1231K"
## [1] "R1232Q"
## [1] "R1232W"
## [1] "K1233E"
## [1] "K1236N"
## [1] "K1236R"
## [1] "V1237F"
## [1] "L1239P"
## [1] "A123G"
## [1] "A123V"
## [1] "E1240Q"
## [1] "D1243N"
## [1] "K1244E"
## [1] "M1245I"
## [1] "T1247I"
## [1] "V1249D"
## [1] "A124D"
## [1] "V1251M"
## [1] "E1253G"
## [1] "L1255M"
## [1] "V125L"
## [1] "A1260D"
## [1] "G1262S"
## [1] "K1264N"
## [1] "K1264R"
## [1] "T1268N"
## [1] "T1268S"
## [1] "N1269S"
## [1] "K126E"
## [1] "A1270S"
## [1] "W1271C"
## [1] "W1273C"
## [1] "D1275N"
## [1] "I1278N"
## [1] "V1279I"
## [1] "V1281F"
## [1] "S1282A"
## [1] "L1283M"
## [1] "A1288G"
## [1] "F1293S"
## [1] "A1294G"
## [1] "E1295K"
## [1] "M1296T"
## [1] "P1298L"
## [1] "R1303Q"
## [1] "R1303W"
## [1] "T1304M"
## [1] "R1306H"
## [1] "R1306S"
## [1] "L1308F"
## [1] "R1309C"
## [1] "R1309H"
## [1] "L1311P"
## [1] "R1316L"
## [1] "R1316Q"
## [1] "F1317C"
## [1] "G1319V"
## [1] "M1320I"
## [1] "R1321K"
## [1] "V1323G"
## [1] "N1325S"
## [1] "A1326S"
## [1] "V1328M"
## [1] "G1329S"
## [1] "A1330D"
## [1] "A1330P"
## [1] "A1330T"
## [1] "I1331V"
## [1] "P1332L"
## [1] "P1332Q"
## [1] "I1334V"
## [1] "M1335R"
## [1] "L1338V"
## [1] "L1339F"
## [1] "p.L1339del"
## [1] "V1340I"
## [1] "F1344L"
## [1] "F1344S"
## [1] "W1345C"
## [1] "L1346I"
## [1] "L1346P"
## [1] "F1348L"
## [1] "N134S"
## [1] "I1350L"
## [1] "I1350T"
## [1] "M1351R"
## [1] "M1351V"
## [1] "V1353M"
## [1] "A1357V"
## [1] "G1358R"
## [1] "G1358W"
## [1] "K1359M"
## [1] "K1359N"
## [1] "M135V"
## [1] "F1360C"
## [1] "R1362S"
## [1] "C1363Y"
## [1] "I1364V"
## [1] "N1365S"
## [1] "Q1366H"
## [1] "Q1366R"
## [1] "G1369R"
## [1] "L136P"
## [1] "D1370A"
## [1] "D1370G"
## [1] "Y1375H"
## [1] "I1377M"
## [1] "V1378M"
## [1] "I137V"
## [1] "N1380K"
## [1] "p.N1380del"
## [1] "S1382I"
## [1] "C1384Y"
## [1] "L1387F"
## [1] "M138I"
## [1] "G1391R"
## [1] "V1398M"
## [1] "p.I137_C139dup"
## [1] "V1400I"
## [1] "V1405L"
## [1] "V1405M"
## [1] "G1406E"
## [1] "G1406R"
## [1] "G1408R"
## [1] "Y1409C"
## [1] "L1412F"
## [1] "Q1414H"
## [1] "A1416E"
## [1] "A1416G"
## [1] "K1419E"
## [1] "I141N"
## [1] "I141V"
## [1] "G1420D"
## [1] "G1420P"
## [1] "G1420R"
## [1] "G1420V"
## [1] "M1422R"
## [1] "D1423H"
## [1] "I1424V"
## [1] "A1427E"
## [1] "A1427S"
## [1] "A1428S"
## [1] "A1428V"
## [1] "D1430N"
## [1] "S1431C"
## [1] "R1432G"
## [1] "R1432S"
## [1] "G1433R"
## [1] "G1433V"
## [1] "G1433W"
## [1] "P1438L"
## [1] "Q1439H"
## [1] "Q1439R"
## [1] "E1441Q"
## [1] "Y1442C"
## [1] "Y1442N"
## [1] "N1443S"
## [1] "L1444I"
## [1] "Y1445H"
## [1] "I1448L"
## [1] "I1448T"
## [1] "Y1449C"
## [1] "Y1449S"
## [1] "V1451D"
## [1] "V1451L"
## [1] "S1458Y"
## [1] "F1460L"
## [1] "T1461S"
## [1] "N1463Y"
## [1] "L1464P"
## [1] "p.F1465_L1480dup"
## [1] "V1468A"
## [1] "V1468F"
## [1] "I1469V"
## [1] "V146A"
## [1] "V146M"
## [1] "N1472S"
## [1] "p.N1472del"
## [1] "F1473C"
## [1] "F1473S"
## [1] "Q1475L"
## [1] "Q1476R"
## [1] "K1477N"
## [1] "K1478E"
## [1] "G1481E"
## [1] "G1481R"
## [1] "G1481V"
## [1] "Q1483H"
## [1] "F1486L"
## [1] "p.F1486del"
## [1] "M1487K"
## [1] "M1487L"
## [1] "T1488R"
## [1] "E1489D"
## [1] "Q1491H"
## [1] "K1493R"
## [1] "p.K1493del"
## [1] "Y1495S"
## [1] "M1498R"
## [1] "M1498T"
## [1] "M1498V"
## [1] "R14C"
## [1] "R14H"
## [1] "R14S"
## [1] "p.K1500del"
## [1] "L1501V"
## [1] "p.L1501_K1505del"
## [1] "G1502A"
## [1] "G1502S"
## [1] "S1503Y"
## [1] "K1504E"
## [1] "K1505N"
## [1] "p.K1505_Q1507del"
## [1] "P1506S"
## [1] "P1506T"
## [1] "p.Q1507_P1509del"
## [1] "P1509T"
## [1] "R1512L"
## [1] "R1512Q"
## [1] "R1512W"
## [1] "L1514M"
## [1] "N1515S"
## [1] "I1521K"
## [1] "I1521T"
## [1] "D1523N"
## [1] "I1524T"
## [1] "V1525A"
## [1] "V1525M"
## [1] "T1526P"
## [1] "K1527R"
## [1] "D152N"
## [1] "V1532F"
## [1] "V1532I"
## [1] "T1533I"
## [1] "C1539F"
## [1] "C1539Y"
## [1] "V1543A"
## [1] "V1543L"
## [1] "T1544P"
## [1] "M1546T"
## [1] "V1547L"
## [1] "E1548K"
## [1] "G1548K"
## [1] "P154L"
## [1] "D1551N"
## [1] "D1551Y"
## [1] "Q1552L"
## [1] "Q1552R"
## [1] "S1553R"
## [1] "E1555K"
## [1] "I1557V"
## [1] "I1559V"
## [1] "L1560F"
## [1] "L1565M"
## [1] "F1567L"
## [1] "A1569P"
## [1] "W156R"
## [1] "I1570V"
## [1] "p.1570_F1571insI"
## [1] "p.I1570dup"
## [1] "F1571C"
## [1] "E1574K"
## [1] "C1575S"
## [1] "T157I"
## [1] "A1581S"
## [1] "L1582P"
## [1] "R1583C"
## [1] "R1583H"
## [1] "Y1585C"
## [1] "F1587V"
## [1] "T1588I"
## [1] "K158T"
## [1] "I1593M"
## [1] "F1594S"
## [1] "F1596C"
## [1] "F1596I"
## [1] "V1597M"
## [1] "V1598A"
## [1] "Y159C"
## [1] "R15G"
## [1] "R15M"
## [1] "R15T"
## [1] "L1601H"
## [1] "I1603F"
## [1] "V1604M"
## [1] "G1605C"
## [1] "G1605D"
## [1] "T1606I"
## [1] "S1609L"
## [1] "S1609W"
## [1] "D1610G"
## [1] "I1611V"
## [1] "Q1613H"
## [1] "Q1613L"
## [1] "p.F1617del"
## [1] "P1619L"
## [1] "P1619Q"
## [1] "E161K"
## [1] "E161Q"
## [1] "T1620K"
## [1] "T1620M"
## [1] "R1623L"
## [1] "R1623Q"
## [1] "V1624I"
## [1] "R1626C"
## [1] "R1626H"
## [1] "R1626L"
## [1] "R1626P"
## [1] "R1629G"
## [1] "R1629Q"
## [1] "Y162C"
## [1] "Y162H"
## [1] "I1630R"
## [1] "I1630V"
## [1] "G1631D"
## [1] "R1632C"
## [1] "R1632H"
## [1] "R1632L"
## [1] "L1634P"
## [1] "R1638Q"
## [1] "G1639A"
## [1] "G1642E"
## [1] "I1643L"
## [1] "R1644C"
## [1] "R1644H"
## [1] "R1644L"
## [1] "T1645M"
## [1] "A1649V"
## [1] "F164L"
## [1] "L1650F"
## [1] "M1652R"
## [1] "M1652T"
## [1] "I1660S"
## [1] "I1660V"
## [1] "G1661E"
## [1] "G1661R"
## [1] "V1667I"
## [1] "M1668T"
## [1] "A166T"
## [1] "S1672Y"
## [1] "F1674V"
## [1] "M1676I"
## [1] "M1676T"
## [1] "N1678S"
## [1] "A1680P"
## [1] "A1680T"
## [1] "Y1681F"
## [1] "W1684R"
## [1] "D1689N"
## [1] "D1690N"
## [1] "A1698T"
## [1] "F16L"
## [1] "M1701I"
## [1] "L1704H"
## [1] "Q1706H"
## [1] "T1708I"
## [1] "p.T1709del"
## [1] "T1709M"
## [1] "T1709R"
## [1] "F170I"
## [1] "S1710L"
## [1] "G1712C"
## [1] "G1712S"
## [1] "D1714G"
## [1] "L1717P"
## [1] "S1718R"
## [1] "N1722D"
## [1] "T1723N"
## [1] "P1725L"
## [1] "C1728R"
## [1] "C1728W"
## [1] "C1728Y"
## [1] "D1729N"
## [1] "P1730A"
## [1] "P1730H"
## [1] "P1730L"
## [1] "G1737D"
## [1] "S1738F"
## [1] "S1738T"
## [1] "R1739Q"
## [1] "R1739W"
## [1] "G1740R"
## [1] "D1741E"
## [1] "D1741N"
## [1] "D1741Y"
## [1] "G1743E"
## [1] "G1743R"
## [1] "S1744I"
## [1] "A1746T"
## [1] "A1746V"
## [1] "V1747M"
## [1] "G1748D"
## [1] "p.G1748del"
## [1] "I1749N"
## [1] "V174I"
## [1] "L1750F"
## [1] "T1753A"
## [1] "I1756V"
## [1] "I1758V"
## [1] "p.I1758del"
## [1] "S1759C"
## [1] "K175N"
## [1] "L1761F"
## [1] "L1761H"
## [1] "I1762M"
## [1] "p.I1762del"
## [1] "V1763L"
## [1] "V1763M"
## [1] "V1764F"
## [1] "M1766L"
## [1] "M1766T"
## [1] "M1766V"
## [1] "Y1767C"
## [1] "I1768V"
## [1] "I1770V"
## [1] "I1771T"
## [1] "L1772V"
## [1] "N1774D"
## [1] "F1775V"
## [1] "V1777L"
## [1] "V1777M"
## [1] "T1779M"
## [1] "L177P"
## [1] "E1780G"
## [1] "E1781D"
## [1] "E1781G"
## [1] "E1784K"
## [1] "L1786Q"
## [1] "L1786R"
## [1] "S1787N"
## [1] "A178G"
## [1] "D1790G"
## [1] "D1790N"
## [1] "p.D1790del"
## [1] "D1792N"
## [1] "D1792V"
## [1] "D1792Y"
## [1] "M1793K"
## [1] "p.Y1795_E1796insD"
## [1] "Y1795C"
## [1] "Y1795H"
## [1] "Y1795N"
## [1] "I1797V"
## [1] "R179Q"
## [1] "I1809M"
## [1] "G180V"
## [1] "Y1811N"
## [1] "S1812L"
## [1] "D1816E"
## [1] "D1816N"
## [1] "D1819N"
## [1] "A1820T"
## [1] "A1820V"
## [1] "E1823K"
## [1] "P1824A"
## [1] "L1825P"
## [1] "R1826C"
## [1] "R1826H"
## [1] "A1828S"
## [1] "A1828T"
## [1] "C182R"
## [1] "C182Y"
## [1] "Q1832E"
## [1] "I1833M"
## [1] "S1834R"
## [1] "L1835F"
## [1] "I1836T"
## [1] "D1839G"
## [1] "M1842L"
## [1] "M1842T"
## [1] "M1842V"
## [1] "G1845R"
## [1] "R1847C"
## [1] "R1847H"
## [1] "H1849R"
## [1] "H184R"
## [1] "C1850S"
## [1] "M1851I"
## [1] "M1851V"
## [1] "D1852V"
## [1] "I1853V"
## [1] "A185T"
## [1] "A185V"
## [1] "V1861F"
## [1] "V1861I"
## [1] "A1870T"
## [1] "K1872N"
## [1] "I1873V"
## [1] "M1875K"
## [1] "M1875T"
## [1] "p.M1875dup"
## [1] "E1877K"
## [1] "T187I"
## [1] "T187S"
## [1] "M1880V"
## [1] "P1884L"
## [1] "Y1889C"
## [1] "E1890K"
## [1] "T1895I"
## [1] "L1896P"
## [1] "R1897Q"
## [1] "R1897W"
## [1] "R1897Y"
## [1] "R1898C"
## [1] "R1898H"
## [1] "R18Q"
## [1] "R18W"
## [1] "E1901K"
## [1] "E1901Q"
## [1] "E1902A"
## [1] "V1903M"
## [1] "S1904L"
## [1] "M1906T"
## [1] "M1906V"
## [1] "I1908V"
## [1] "Q1909R"
## [1] "R190G"
## [1] "R190L"
## [1] "R190Q"
## [1] "R190W"
## [1] "R1910K"
## [1] "R1913C"
## [1] "R1913H"
## [1] "R1913S"
## [1] "R1914G"
## [1] "H1915N"
## [1] "H1915P"
## [1] "H1915Q"
## [1] "H1915Y"
## [1] "R1919C"
## [1] "R1919H"
## [1] "S1920C"
## [1] "K1922N"
## [1] "K1922R"
## [1] "H1923D"
## [1] "H1923Y"
## [1] "A1924T"
## [1] "S1925F"
## [1] "L1927P"
## [1] "F1928V"
## [1] "R1929C"
## [1] "R1929H"
## [1] "Q1930H"
## [1] "A1932V"
## [1] "G1933A"
## [1] "G1933D"
## [1] "G1933V"
## [1] "G1935S"
## [1] "S1937A"
## [1] "E1938K"
## [1] "p.E1939_E1943del"
## [1] "W193R"
## [1] "P1942H"
## [1] "P1942S"
## [1] "R1944Q"
## [1] "I1948V"
## [1] "A1949S"
## [1] "A1949T"
## [1] "Y1950C"
## [1] "V1951L"
## [1] "V1951M"
## [1] "E1954K"
## [1] "N1955Y"
## [1] "S1957P"
## [1] "R1958P"
## [1] "R1958Q"
## [1] "P1962L"
## [1] "P1962S"
## [1] "P1963L"
## [1] "S1964F"
## [1] "S1965G"
## [1] "S1965N"
## [1] "I1968M"
## [1] "I1968N"
## [1] "I1968S"
## [1] "I1968T"
## [1] "I1968V"
## [1] "p.S1970_S1972del"
## [1] "F1973L"
## [1] "P1975T"
## [1] "S1976C"
## [1] "Y1977N"
## [1] "V1980F"
## [1] "R1982T"
## [1] "A1983G"
## [1] "A1983V"
## [1] "T1984I"
## [1] "S1985R"
## [1] "D1986G"
## [1] "D1986N"
## [1] "N1987K"
## [1] "L1988R"
## [1] "V1990L"
## [1] "R1991Q"
## [1] "R1991W"
## [1] "G1992A"
## [1] "S1996N"
## [1] "S1996R"
## [1] "H1997R"
## [1] "S199T"
## [1] "M1I"
## [1] "D2000Y"
## [1] "A2002T"
## [1] "D2003N"
## [1] "F2004I"
## [1] "F2004L"
## [1] "F2004V"
## [1] "p.F2004dup"
## [1] "P2005A"
## [1] "P2005L"
## [1] "P2005S"
## [1] "p.Pro2006del"
## [1] "P2006A"
## [1] "P2006T"
## [1] "P2008L"
## [1] "D2009E"
## [1] "R2010G"
## [1] "R2012C"
## [1] "R2012H"
## [1] "V2016L"
## [1] "V2016M"
## [1] "I202T"
## [1] "A204T"
## [1] "A204V"
## [1] "E208K"
## [1] "N209S"
## [1] "N209T"
## [1] "S20F"
## [1] "I210T"
## [1] "L212P"
## [1] "L212Q"
## [1] "S216L"
## [1] "R219C"
## [1] "R219H"
## [1] "L21V"
## [1] "T220I"
## [1] "R222L"
## [1] "R222Q"
## [1] "V223L"
## [1] "L224F"
## [1] "R225Q"
## [1] "R225W"
## [1] "A226G"
## [1] "A226V"
## [1] "L227P"
## [1] "K228R"
## [1] "A22V"
## [1] "I230M"
## [1] "I230T"
## [1] "I230V"
## [1] "V232F"
## [1] "V232I"
## [1] "P234S"
## [1] "G235R"
## [1] "I239V"
## [1] "I239V "
## [1] "A23S"
## [1] "V240M"
## [1] "A242D"
## [1] "Q245K"
## [1] "V247L"
## [1] "V258A"
## [1] "E25K"
## [1] "S262G"
## [1] "V263I"
## [1] "A265V"
## [1] "L266H"
## [1] "G268S"
## [1] "Q270K"
## [1] "L271V"
## [1] "G274C"
## [1] "N275K"
## [1] "L276P"
## [1] "L276Q"
## [1] "H278D"
## [1] "H278R"
## [1] "R27C"
## [1] "R27H"
## [1] "R27L"
## [1] "C280Y"
## [1] "V281M"
## [1] "R282C"
## [1] "R282H"
## [1] "T285K"
## [1] "A286S"
## [1] "A286V"
## [1] "N288S"
## [1] "G289S"
## [1] "M28I"
## [1] "M28L"
## [1] "M28T"
## [1] "p.T290_G292del"
## [1] "p.N291_S293dup"
## [1] "N291H"
## [1] "N291K"
## [1] "N291S"
## [1] "G292S"
## [1] "V294M"
## [1] "E295K"
## [1] "G298D"
## [1] "G298S"
## [1] "L299F"
## [1] "L299M"
## [1] "A29E"
## [1] "A29V"
## [1] "A2T"
## [1] "V300I"
## [1] "D305N"
## [1] "L306F"
## [1] "L306V"
## [1] "S309C"
## [1] "S309N"
## [1] "E30G"
## [1] "D310N"
## [1] "E312K"
## [1] "L315P"
## [1] "K317E"
## [1] "K317M"
## [1] "K317N"
## [1] "G319C"
## [1] "G319R"
## [1] "G319S"
## [1] "T320N"
## [1] "S321Y"
## [1] "L325R"
## [1] "S330F"
## [1] "A332T"
## [1] "C335R"
## [1] "C335S"
## [1] "P336L"
## [1] "R340Q"
## [1] "R340W"
## [1] "C341Y"
## [1] "A344S"
## [1] "E346D"
## [1] "E346G"
## [1] "E346K"
## [1] "P348A"
## [1] "D349N"
## [1] "R34C"
## [1] "R34H"
## [1] "H350Q"
## [1] "G351C"
## [1] "G351D"
## [1] "G351S"
## [1] "G351V"
## [1] "Y352C"
## [1] "T353I"
## [1] "F355C"
## [1] "F355I"
## [1] "D356N"
## [1] "A359T"
## [1] "G35S"
## [1] "R367C"
## [1] "R367H"
## [1] "R367L"
## [1] "M369K"
## [1] "T370M"
## [1] "Q371E"
## [1] "W374G"
## [1] "R376C"
## [1] "R376H"
## [1] "T37A"
## [1] "S384T"
## [1] "A385T"
## [1] "G386E"
## [1] "G386R"
## [1] "I388S"
## [1] "Y389H"
## [1] "V396A"
## [1] "V396L"
## [1] "I397F"
## [1] "I397T"
## [1] "I397V"
## [1] "L39F"
## [1] "N3K"
## [1] "N3S"
## [1] "G400A"
## [1] "G400E"
## [1] "G400R"
## [1] "S401P"
## [1] "F402L"
## [1] "L404Q"
## [1] "L404V"
## [1] "N406K"
## [1] "N406S"
## [1] "L409P"
## [1] "L409V"
## [1] "A410V"
## [1] "V411M"
## [1] "V412D"
## [1] "A413E"
## [1] "A413T"
## [1] "M414V"
## [1] "A415T"
## [1] "Y416C"
## [1] "E418K"
## [1] "I424M"
## [1] "A425P"
## [1] "A425T"
## [1] "E428K"
## [1] "E429K"
## [1] "p.E429del"
## [1] "K430E"
## [1] "R433C"
## [1] "R433H"
## [1] "R433S"
## [1] "A437V"
## [1] "M438L"
## [1] "M438T"
## [1] "E439K"
## [1] "E439V"
## [1] "R43Q"
## [1] "L441F"
## [1] "H445D"
## [1] "H445Q"
## [1] "H445Y"
## [1] "E446K"
## [1] "A447G"
## [1] "A447S"
## [1] "T449A"
## [1] "Y449C"
## [1] "G452D"
## [1] "V453M"
## [1] "V456M"
## [1] "R458C"
## [1] "R458H"
## [1] "S459G"
## [1] "G45A"
## [1] "L461V"
## [1] "E462A"
## [1] "E462K"
## [1] "M463R"
## [1] "M463T"
## [1] "L466F"
## [1] "P468L"
## [1] "V469I"
## [1] "N470K"
## [1] "R474G"
## [1] "R474K"
## [1] "R475K"
## [1] "R475S"
## [1] "K480N"
## [1] "R481Q"
## [1] "R481W"
## [1] "M482I"
## [1] "T486A"
## [1] "T486S"
## [1] "E48K"
## [1] "G490A"
## [1] "G490E"
## [1] "E491G"
## [1] "R493K"
## [1] "K496M"
## [1] "K496N"
## [1] "S497C"
## [1] "E49K"
## [1] "F4V"
## [1] "E500K"
## [1] "D501G"
## [1] "P503S"
## [1] "R504T"
## [1] "A505E"
## [1] "M506K"
## [1] "p.N507_L515dup"
## [1] "T512I"
## [1] "R513C"
## [1] "R513H"
## [1] "R513P"
## [1] "G514C"
## [1] "R517S"
## [1] "S519F"
## [1] "A51V"
## [1] "M520R"
## [1] "M520V"
## [1] "K521E"
## [1] "P522S"
## [1] "R523C"
## [1] "R523H"
## [1] "R523S"
## [1] "S524Y"
## [1] "S525G"
## [1] "R526C"
## [1] "R526H"
## [1] "G527R"
## [1] "S528R"
## [1] "P52H"
## [1] "P52S"
## [1] "F530S"
## [1] "F530V"
## [1] "T531A"
## [1] "T531I"
## [1] "F532C"
## [1] "F532L"
## [1] "R533C"
## [1] "R533H"
## [1] "R533S"
## [1] "R535G"
## [1] "R535Q"
## [1] "D536H"
## [1] "G538D"
## [1] "R53Q"
## [1] "F543L"
## [1] "A551T"
## [1] "A551V"
## [1] "G552R"
## [1] "G552W"
## [1] "E553K"
## [1] "S554I"
## [1] "S554N"
## [1] "E555K"
## [1] "H557L"
## [1] "H557Q"
## [1] "H557Y"
## [1] "H558R"
## [1] "T559I"
## [1] "T559R"
## [1] "V563G"
## [1] "L567Q"
## [1] "R568C"
## [1] "R568H"
## [1] "R569Q"
## [1] "R569W"
## [1] "T570N"
## [1] "S571I"
## [1] "A572D"
## [1] "A572F"
## [1] "A572S"
## [1] "A572V"
## [1] "Q573E"
## [1] "Q573R"
## [1] "G574E"
## [1] "S577N"
## [1] "P578R"
## [1] "P578T"
## [1] "G579R"
## [1] "S581L"
## [1] "P583L"
## [1] "G584R"
## [1] "A586G"
## [1] "A586T"
## [1] "p.A586_L587del"
## [1] "p.586_587delAL"
## [1] "H588N"
## [1] "H588R"
## [1] "K590Q"
## [1] "N592K"
## [1] "N592S"
## [1] "D596G"
## [1] "C597G"
## [1] "C597Y"
## [1] "G599R"
## [1] "V601A"
## [1] "L604V"
## [1] "G605E"
## [1] "A606T"
## [1] "G607D"
## [1] "G607R"
## [1] "G607V"
## [1] "D608H"
## [1] "D608N"
## [1] "A60P"
## [1] "E610K"
## [1] "P614S"
## [1] "G615E"
## [1] "S616N"
## [1] "L618F"
## [1] "L619F"
## [1] "R620C"
## [1] "R620H"
## [1] "M623T"
## [1] "L624Q"
## [1] "E625D"
## [1] "E625Q"
## [1] "P627L"
## [1] "P628R"
## [1] "D629Y"
## [1] "K62T"
## [1] "T630M"
## [1] "T632M"
## [1] "S634L"
## [1] "S634W"
## [1] "E636K"
## [1] "P637L"
## [1] "G638D"
## [1] "G639A"
## [1] "G639R"
## [1] "K63N"
## [1] "P640A"
## [1] "P640L"
## [1] "P640S"
## [1] "A647D"
## [1] "A647S"
## [1] "A647V"
## [1] "P648L"
## [1] "C649R"
## [1] "C649Y"
## [1] "p.L64del"
## [1] "D651H"
## [1] "G652D"
## [1] "G652S"
## [1] "E654D"
## [1] "E654K"
## [1] "E654Q"
## [1] "E655K"
## [1] "P656L"
## [1] "A658V"
## [1] "R659Q"
## [1] "R659W"
## [1] "R661Q"
## [1] "R661W"
## [1] "A662S"
## [1] "S664G"
## [1] "A665S"
## [1] "A665T"
## [1] "V668I"
## [1] "D66N"
## [1] "S671C"
## [1] "A672S"
## [1] "A672T"
## [1] "L673P"
## [1] "L673V"
## [1] "E677V"
## [1] "R680C"
## [1] "R680H"
## [1] "H681P"
## [1] "C683G"
## [1] "C683R"
## [1] "C683S"
## [1] "C686F"
## [1] "R689C"
## [1] "R689H"
## [1] "Y68C"
## [1] "A691S"
## [1] "A691T"
## [1] "Q692K"
## [1] "R693C"
## [1] "R693H"
## [1] "Y694C"
## [1] "I696N"
## [1] "G69D"
## [1] "L6S"
## [1] "P701L"
## [1] "M704T"
## [1] "S705F"
## [1] "S705P"
## [1] "G709V"
## [1] "N70K"
## [1] "N70S"
## [1] "V714A"
## [1] "V714D"
## [1] "M715I"
## [1] "M715T"
## [1] "P717L"
## [1] "F718L"
## [1] "I723V"
## [1] "T724I"
## [1] "V728I"
## [1] "p.L729del"
## [1] "N730K"
## [1] "T731I"
## [1] "F733L"
## [1] "M734V"
## [1] "A735E"
## [1] "A735T"
## [1] "A735V"
## [1] "L736P"
## [1] "p.N740del"
## [1] "p.M741_T742delinsI "
## [1] "T742A"
## [1] "E746K"
## [1] "E747A"
## [1] "M748I"
## [1] "E74D"
## [1] "E74K"
## [1] "Q750R"
## [1] "V751F"
## [1] "V751I"
## [1] "G752R"
## [1] "G758E"
## [1] "I759V"
## [1] "E763D"
## [1] "E763K"
## [1] "M764K"
## [1] "M764R"
## [1] "L771V"
## [1] "D772N"
## [1] "P773S"
## [1] "Y774C"
## [1] "Y774D"
## [1] "p.Y776del"
## [1] "F777L"
## [1] "Q779K"
## [1] "G77R"
## [1] "N782T"
## [1] "I783T"
## [1] "F784L"
## [1] "D785N"
## [1] "I788V"
## [1] "V789A"
## [1] "V789I"
## [1] "L791F"
## [1] "L793F"
## [1] "G797V"
## [1] "P79A"
## [1] "P79H"
## [1] "P79S"
## [1] "P7L"
## [1] "P7R"
## [1] "R800C"
## [1] "R800H"
## [1] "R800L"
## [1] "M801V"
## [1] "p.801_803delMSN/insS"
## [1] "S805L"
## [1] "V806M"
## [1] "R808C"
## [1] "R808H"
## [1] "R808P"
## [1] "R811G"
## [1] "R811H"
## [1] "L812Q"
## [1] "R814Q"
## [1] "F816L"
## [1] "F816Y"
## [1] "K817E"
## [1] "W822C"
## [1] "P823T"
## [1] "N826D"
## [1] "L828V"
## [1] "D82E"
## [1] "G833R"
## [1] "N834D"
## [1] "S835A"
## [1] "S835L"
## [1] "V836M"
## [1] "L839P"
## [1] "N841K"
## [1] "T843A"
## [1] "L846R"
## [1] "I848F"
## [1] "D84G"
## [1] "D84N"
## [1] "V850M"
## [1] "F851L"
## [1] "V856L"
## [1] "G857D"
## [1] "M858L"
## [1] "P85S"
## [1] "P85T"
## [1] "S866L"
## [1] "S866P"
## [1] "E867K"
## [1] "E867Q"
## [1] "R869S"
## [1] "F86L"
## [1] "D872N"
## [1] "S873A"
## [1] "G874D"
## [1] "R878C"
## [1] "R878H"
## [1] "R878L"
## [1] "W879R"
## [1] "Y87C"
## [1] "H886P"
## [1] "H886Q"
## [1] "S88G"
## [1] "I890T"
## [1] "I891N"
## [1] "I891T"
## [1] "F892I"
## [1] "R893C"
## [1] "R893H"
## [1] "I894M"
## [1] "L895F"
## [1] "C896S"
## [1] "G897E"
## [1] "G897R"
## [1] "R8P"
## [1] "R8Q"
## [1] "R8W"
## [1] "E901K"
## [1] "S901L"
## [1] "S910L"
## [1] "G911E"
## [1] "Q912R"
## [1] "C915R"
## [1] "L917R"
## [1] "L917V"
## [1] "V922I"
## [1] "V924I"
## [1] "I925F"
## [1] "N927K"
## [1] "N927S"
## [1] "L928P"
## [1] "T92I"
## [1] "L935P"
## [1] "L939F"
## [1] "F93S"
## [1] "S940N"
## [1] "S941F"
## [1] "S941N"
## [1] "S943N"
## [1] "D945G"
## [1] "I94S"
## [1] "I94V"
## [1] "E952K"
## [1] "D953E"
## [1] "D953Y"
## [1] "R954G"
## [1] "M956I"
## [1] "M956T"
## [1] "L959P"
## [1] "V95I"
## [1] "V95L"
## [1] "Q960K"
## [1] "A964G"
## [1] "R965C"
## [1] "R965H"
## [1] "R965L"
## [1] "Q967R"
## [1] "G969C"
## [1] "G969S"
## [1] "R971C"
## [1] "R971H"
## [1] "K974D"
## [1] "R975Q"
## [1] "R975W"
## [1] "D979H"
## [1] "C981F"
## [1] "C982R"
## [1] "G983D"
## [1] "R986L"
## [1] "R986Q"
## [1] "R986W"
## [1] "R988Q"
## [1] "R988W"
## [1] "K991E"
## [1] "K991T"
## [1] "A993S"
## [1] "A993T"
## [1] "A994T"
## [1] "A994V"
## [1] "L995F"
## [1] "A996T"
## [1] "A997D"
## [1] "A997S"
## [1] "A997T"
## [1] "G999D"
## [1] "G9S"
## [1] "G9V"

Back to top

Step 6 Fit posterior mean \(\tilde\mu_i\) with clinical covariates

The clinical covariates include eaRate, blastpssm, provean_score, pph2_prob, ipeak, and penetrance density.

Here pattern mixture model is used.

Back to top

Step 7 Use fitted penetrance and MSE from Step 6 to solve for \(\alpha\) and \(\beta\)

Adjust \(\alpha_{fi}\) and \(\beta_{fi}\) to \(\alpha_0\) and \(\beta_0\) if they were too small (<0.01)

Back to top

Step 8 Calculate convergence criterion \(\delta\) to see if it’s small enough to end the loop

Calculate the convergence criterion \(\delta\) as follows, and if it’s greater than 0.1, then repeat Step 6 to Step 8. When \(\delta\) is smaller than 0.1, EM converges.

\[\delta=\sum_{i=1}^m |p^{em+}_i-p^{em-}_i| \]

where \(p_i^{em+}\) is the penetrance after one iteration of EM alogrithm for each variant, and \(p_i^{em-}\) is the penetrance from the previous iteration of EM alogrithm for each variant.

EM algorithm converges in 7 steps.

delta <- 10

d$eaRate <- as.numeric(d$eaRate)
d$blastpssm <- as.numeric(d$blastpssm)
d$ipeak <- suppressWarnings(as.numeric(d$ipeak))

covariates <- c("eaRate","blastpssm","provean_score","pph2_prob","ipeak","feat_dist_w")

# weight for weighted regression
d$weight_r <- 1/p_variance_w

d$alpha_1_w <- alpha_1_w
d$beta_1_w <- beta_1_w
d$p_mean_w <- p_mean_w

regression <- function(dv, pivs, nivs, data) {
  # run a linear model with text arguments for dv and ivs
  piv_string <- paste(pivs, collapse=" + ")
  niv_string <- paste(nivs, collapse=" - ")
  if(niv_string!="") iv_string <- paste(piv_string, " - ", niv_string, sep = "")
  if(niv_string=="") iv_string <- paste(piv_string)
  #print(iv_string)
  regression_formula <- as.formula(paste(dv, iv_string, sep=" ~ "))
  #print(regression_formula)
  lm(regression_formula, data, weights = d$weight_r)
}

count <- 0

while(delta>0.1 & count < 10){
  
  count <- count + 1
  alpha_f <- NULL
  beta_f <- NULL
  
  for(i in 1:nrow(d)){
  newdata = data.frame(var=d[i,"var"])
  newdata[covariates] <- d[i,covariates]
  model <- regression("p_mean_w", covariates, 
                      colnames(newdata)[colSums(is.na(newdata))>0], d)
  mean_f <- predict(model, newdata)
  variance_f <- (predict(model, newdata,se.fit = T)$se.fit)^2
  alpha <- solab(mean_f,variance_f)[1]
  beta <- solab(mean_f,variance_f)[2]
  if(alpha<0.01 | beta<0.01){
    alpha_f[i]=alpha_0_w
    beta_f[i]=beta_0_w
  }else{
    alpha_f[i]=alpha
    beta_f[i]=beta
  }
  }
  
  new_mean <- (alpha_f + d$brs1)/(alpha_f + beta_f + d$total_carriers)
  
  delta <- sum(abs(new_mean-d$p_mean_w))
  
  d$p_mean_w <- new_mean

}

Back to top

Scale the variance by a factor of 20

Results

Part 1 Coverage plots

Bootstrap and get the coverage rate

  1. Use the observed penetrance from some variant as the TRUE penetrance for that variant, generate n binomial observations

  2. Use the final EM algorithm posterior as the prior for Beta-Binomial, incorporate data from previous step, generate the posterior distribution, get 95% credible interval.

  3. Check whether the interval cover the true penetrance from Step 1.

  4. Repeat Step 1 to Step 3 N times to get the coverage rate.

Plot coverage

Part 2 Output Bland-Altman plots

Part 4 Correlation Coefficients and Plots

The code to produce the weighted correlations is shown in the code chunk. The size of the circles in the presented figures are scaled by the Log10(Total Number of Carriers).

## wCorr v1.9.1

## [1] "Peak Current Covariate Alone"
## [1] 0.3519097
## [1] 0.3519097
##      2.5%     97.5% 
## 0.2436698 0.4620624
## [1]  2.00000 17.79004

## [1] "Penetrance Density Covariate Alone"
## [1] 0.6617035
## [1] 0.6617035
##      2.5%     97.5% 
## 0.5228591 0.7578649
## [1]    2.0000 -128.4841

## [1] "Peak Current and Penetrance Density Covariates"
## [1] 0.7562233
## [1] 0.7562233
##      2.5%     97.5% 
## 0.6540318 0.8305172
## [1]    3.0000 -200.2098

## [1] "In Silico, Peak Current, and Penetrance Density Covariates"
## [1] 0.7806726
## [1] 0.7806726
##      2.5%     97.5% 
## 0.6937503 0.8544453
## [1]    7.0000 -215.9894

## [1] "In Silico Covariates Alone"
## [1] 0.1942237
## [1] 0.1942237
##      2.5%     97.5% 
## 0.1244528 0.2776599
## [1]  5.00000 72.78969

## [1] "Peak Current Covariate Alone"
## [1] 0.2226512
## [1] 0.2226512
##      2.5%     97.5% 
## 0.1296066 0.3266152
## [1]   2.0000 155.3884

## [1] "Penetrance Density Covariate Alone"
## [1] 0.3513447
## [1] 0.3513447
##      2.5%     97.5% 
## 0.1826644 0.4897752
## [1]   2.0000 114.6661

## [1] "Peak Current and Penetrance Density Covariates"
## [1] 0.4225387
## [1] 0.4225387
##      2.5%     97.5% 
## 0.2794866 0.5603973
## [1]  3.00000 90.50763

## [1] "In Silico, Peak Current, and Penetrance Density Covariates"
## [1] 0.436451
## [1] 0.436451
##      2.5%     97.5% 
## 0.2783023 0.5644973
## [1]  7.00000 93.02051

## [1] "In Silico Covariates Alone"
## [1] 0.1174166
## [1] 0.1174166
##       2.5%      97.5% 
## 0.05564675 0.19489625
## [1]   5.0000 189.9553